Неделя один

Обозначения

Одна зависимая (объясняемая) переменная: y
Несколько регрессоров (предикторов, объясняющих переменных): x, z…
По каждой переменной n наблюдений: $y_1, y_2...y_n$

Модель — это формуля для объясняемой переменной.

Пример

Возьмём например данные по машинам 1920-х годов. Тут видимая линейная взаимосвязь.

Данные по машинам 1920-х гдов

Модель для этих данных может иметь вид $y_i = \beta_1 + \beta_2x_i + \varepsilon_i$

Что здесь есть что? У нас есть наблюдаемые переменные,

y — это длина тормозного пути
x — это скорость, с которой ехала машина.
Есть неизвестные коэффициенты $ _1, _2$
случайная составляющая, ошибка

То есть, $\beta_2$ показывает — насколько увеличивается тормозной путь, если машина разгонится на один лишний километр в час. И есть некая случайная составляющая $ε_i$, это может быть все что угодно:

водитель по-другому нажимал на тормоз,
что-то там на дороге было другое,

то есть это та часть, по которой у нас нет возможности предсказать, но, тем не менее, вот эта случайная ошибка она входит в y.

План действий

Придумать адекватную модель
Получить оценки неизвестных параметров $\widehat\beta_1, \widehat\beta_2$
Спрогнозировать, заменив неизвестные параметры на оценки $\widehat y_i = \widehat\beta_1 + \widehat\beta_2 x_i$

Метод наименьших квадратов

Как найти $\beta_1, \beta_2$? Собственно методом наименьших квадратов.

Если мы придумали какие-то оценки, $\beta_1, \beta_2$ то соответственно возникает такое понятие как ошибка прогноза

\[\widehat \varepsilon_i = y_i - \widehat y_i\]

Есть суммарная ошибка, чтобы суммарная ошибка не занулялась одна в плюс, другая в минус, не компенсировали друг друга, мы возведем в квадрат. И посчитаем сумму квадратов ошибок прогноза, то есть сумма $\widehat \varepsilon_i ^2$

\[Q(\widehat\beta_1, \widehat\beta_2) = \sum_{i=1}^{n} \widehat\varepsilon^{2}_{i} = \sum_{i=1}^{n} (y_i - \widehat y_i)^2 \]

Суть метода МНК

Возьмите в качестве оценок такие коэффициенты $\widehat\beta_1, \widehat\beta_2$ при которых сумма квадратов ошибок прогноза будет минимальна

Случай для одного регрессора

Всё сводится к тому, что в модели \[M_1: y_i = \beta + \varepsilon_i\]

\[\widehat\beta = \bar y\]

Случай для одного регрессора

Случай для двух регрессоров

Тут несколько сложнее, подробности на скриншоте

Случай для двух регрессоров

Случай для множества регрессоров

Что мы получили:

Случай для двух регрессоров

Ещё раз пробежимся по терминам

Случай для двух регрессоров продолжение

Графическое представление

Изобразим на графике, всё то, о чем мы только что проговорили.

Графический вывод

В частности $\widehat\beta_1$, т.е. точка на оси ординат от пересечении с прямой регрессии, называют пересечением или Intercept-ом

Метод наименьших квадратов подбирает прямую так, чтобы суммарные расстояния от точек до прямой были минимальными.

Графический вывод продолжение

Множественные регрессоры

Случай множественных регрессоров принципиально не отличается от двух регрессоров. Поэтому рассмотрим на примере трёх предикторов.

Множественные регрессоры

Суммы квадратов

Что есть что:

RSS (Residual Sum of Squares) — этот показатель меряет, насколько велики эпсилон с крышкой, насколько они далеко лежат от нуля.
TSS (Total Sum of Squares) — этот показатель меряет, насколько каждый из $y_i$ не похож на $y$ среднее. Если $y_i$ далеко лежит от $y$ среднее, соответственно, это слагаемое в сумме квадратов будет большим. И вся сумма квадратов будет большой.
ESS (Explained Sum of Squares) — она показывает, насколько прогнозное значение $y_i$ с крышкой далеко легло от $y$ среднего.

Лекбез по линейной алгебре

Язык эконометрики во многом — это язык линейной алгебры, нужно его знать.

Обозначения

Маленькой буквой $y$ мы будем обозначать вектор, то есть столбик из всех игреков, записанных друг под другом: у_1, у_2 и так далее, у_n.
Ну, соответственно, $х$ маленькое — это все иксы: х_1, х_2 и так далее, х_n, записанные друг под другом.
Аналогично $\widehat\beta$.
И еще введем обозначение: единичку со стрелочкой. Это будет вектор-столбец, то есть столбик из единичек в количестве n штук, потому что у нас в модели будет n наблюдений.

Тогда для нашей модели

\[\widehat y = \widehat\beta_1 \cdot \overrightarrow 1 + \widehat\beta_2 \cdot x + \widehat\beta_3 \cdot z\]

Большая буква X — это все регрессоры, помещенные в одну большую табличку чисел, которые называются матрицей.

Таблица регрессоров

Рассмотрим ещё такое понятие как длина вектора

Длина вектора

Где у нас бывают длины векторов:

Пример длины вектора

Есть одно замечательное применение. Если скалярное произведение векторов равно нулю, значит они перпендикулярны.

Скалярное произведение векторов

x <- c(1,2, -3)
y <- c(-6,0,-2)

sum(x * y)

## [1] 0

Линейная модель регрессии в n-мерном пространстве

n-мерное простарнство

Есть вектор y, есть вектор из одних 1, я продолжаю этот вектор до прямой и оказывается, что

в простой модели без регрессоров спроецировав вектор y на эту прямую, я получу вектор y с крышкой — предсказанные значения.

Соответственно, мы получили попутно еще один факт: если любой вектор проецировать на вектор из 1, то получится вектор средних значений

Геометрия множественной регрессии

Напомню, что мы вывели условие первого порядка

Геометрическая интерпретация условий первого порядка

Это означает, что вектора перпендикулярны.

Облачко это все те вектора, которые можно получить, складывая с некоторыми весами вектор x, вектор z и вектор из единичек. Это гиперплоскость.

Мы получаем следующую интрпретацию метода наименьших квадратов:

Первый геометрический факт Прогнозы, которые мы получаем по методу наименьших квадратов — это проекция исходного вектора зависимых переменных y на множество векторов, получаемых с помощью сложения с разными весами вектора из единичек, вектора x и вектора z.
Второй геометрический факт Если я спроецировать вектор y на прямую, порождаемую вектором из единичек, и спроецировать y с крышечкой на ту же самую прямую, то эти проекции попадут в одну и ту же точку.
$TSS = ESS + RSS$
$\frac{ESS}{TSS}= \frac{BC^2}{AB^2}=(\frac{BC}{AB})^2= (cos \varphi)^2$

Геометрчиеские факты

Коэффициент детерминации

Если в регрессию включён свободный член $\beta_1$ то действуют следующие правила:

Правила при включении свободного члена в регрессию

Эти правила позволяют придумать простой показатель качества работы модели — коэффициент детерминации.

Чем прогнозы точнее похожи на настоящие значения, тем меньше будут ошибки прогнозов и тем меньше будет сумма квадратов ошибок прогнозов RSS. Соответственно $\frac{ESS}{TSS} = R^2$ будет примерно равен 1 если RSS будет у ноля. Соответственно мы получим коэффициент детерминации, который всегда лежит от 0 до 1.

Другими словами — коэффициент детерминации это доля объяснённого разброса в общем разбросе.

Правила при включении свободного члена в регрессию

С одной стороны, коэффициент детерминации — это доля объясненной дисперсии игрек, доля объясненного разброса.

С другой стороны, коэффициент детерминации — это выборочная корреляция между прогнозами и настоящим игрек, взятое в квадрат.

Чем коэффициент детерминации выше, тем больше предсказание похожи на реальные значения. Чем коэффициент детерминации выше, тем выше доля объясненной дисперсии.

Пример с фертильностью. МНК в R

Посмотрим как зависит фертильность от других

t <- swiss
# нарисовать много диаграм рассеивания
# install.packages("GGally") # матрица диаграмм рассеяния
library("GGally")

## Loading required package: ggplot2

str(t)

## 'data.frame':    47 obs. of  6 variables:
##  $ Fertility       : num  80.2 83.1 92.5 85.8 76.9 76.1 83.8 92.4 82.4 82.9 ...
##  $ Agriculture     : num  17 45.1 39.7 36.5 43.5 35.3 70.2 67.8 53.3 45.2 ...
##  $ Examination     : int  15 6 5 12 17 9 16 14 12 16 ...
##  $ Education       : int  12 9 5 7 15 7 7 8 7 13 ...
##  $ Catholic        : num  9.96 84.84 93.4 33.77 5.16 ...
##  $ Infant.Mortality: num  22.2 22.2 20.2 20.3 20.6 26.6 23.6 24.9 21 24.4 ...

ggpairs(t)

Уже тут можно посмотреть много всяких корреляций.

Перейдём к оценке

model2 <-  lm(data = t, Fertility ~ Agriculture + Education + Catholic)
# КОэффициенты
coef(model2)

## (Intercept) Agriculture   Education    Catholic 
##  86.2250198  -0.2030377  -1.0721468   0.1452013

# Спрогнозированные значения
fitted(model2)

##   Courtelary     Delemont Franches-Mnt      Moutier   Neuveville 
##     71.35382     79.73758     86.36550     76.21257     62.05992 
##   Porrentruy        Broye        Glane      Gruyere       Sarine 
##     84.70365     77.94869     77.98965     82.07990     76.37831 
##      Veveyse        Aigle      Aubonne     Avenches     Cossonay 
##     81.01451     62.00804     65.34456     61.67811     67.20324 
##    Echallens     Grandson     Lausanne    La Vallee       Lavaux 
##     72.85406     71.22373     54.02437     62.00809     62.16632 
##       Morges       Moudon        Nyone         Orbe         Oron 
##     64.12130     72.47751     65.22299     69.41765     71.04507 
##      Payerne Paysd'enhaut        Rolle        Vevey      Yverdon 
##     66.61076     70.48740     64.27982     63.09324     68.48321 
##      Conthey    Entremont       Herens     Martigwy      Monthey 
##     81.11782     77.02791     80.38838     78.28372     84.09311 
##   St Maurice       Sierre         Sion       Boudry La Chauxdfnd 
##     75.54878     80.27332     73.53528     66.37864     74.87034 
##     Le Locle    Neuchatel   Val de Ruz ValdeTravers V. De Geneve 
##     70.52554     50.79966     71.80743     76.17918     35.30542 
##  Rive Droite  Rive Gauche 
##     52.99371     57.97821

# Остатки
residuals(model2)

##   Courtelary     Delemont Franches-Mnt      Moutier   Neuveville 
##    8.8461774    3.3624191    6.1345049    9.5874340   14.8400828 
##   Porrentruy        Broye        Glane      Gruyere       Sarine 
##   -8.6036475    5.8513084   14.4103471    0.3201012    6.5216935 
##      Veveyse        Aigle      Aubonne     Avenches     Cossonay 
##    6.0854871    2.0919628    1.5554442    7.2218873   -5.5032423 
##    Echallens     Grandson     Lausanne    La Vallee       Lavaux 
##   -4.5540632    0.4762715    1.6756342   -7.7080933    2.9336804 
##       Morges       Moudon        Nyone         Orbe         Oron 
##    1.3786987   -7.4775133   -8.6229883  -12.0176461    1.4549265 
##      Payerne Paysd'enhaut        Rolle        Vevey      Yverdon 
##    7.5892409    1.5125978   -3.7798150   -4.7932370   -3.0832083 
##      Conthey    Entremont       Herens     Martigwy      Monthey 
##   -5.6178154   -7.7279097   -3.0883806   -7.7837172   -4.6931098 
##   St Maurice       Sierre         Sion       Boudry La Chauxdfnd 
##  -10.5487835   11.9266828    5.7647206    4.0213575   -9.1703410 
##     Le Locle    Neuchatel   Val de Ruz ValdeTravers V. De Geneve 
##    2.1744592   13.6003353    5.7925740   -8.5791790   -0.3054172 
##  Rive Droite  Rive Gauche 
##   -8.2937095  -15.1782122

# ПОказатель RSS
deviance(model2)

## [1] 2567.884

# Показатель R квадрат
report <- summary(model2)
report$r.squared

## [1] 0.6422541

# Это тоже самое что
cor(t$Fertility, fitted(model2))^2

## [1] 0.6422541

Неделя два

Вторая неделя — стат.свойства оценок коэффициентов

Формулируем стандартные предпосылки — они ведут к свойствам — свойства позволят строить доверительные интервалы.

Для формулировки предпосылок — сформулируем поняти условное мат.ожидание.

Условное математическое ожидание

Условия мат. ожидания

Если формально: это некая случайная величина s с тильдой, которая является, с одной стороны, функцией от r, и, с другой стороны, эта самая s с тильдой очень похожа на s, а именно: математической ожидание от s с тильдой равно математическому ожиданию от s, и ковариация между s и любой функцией от r равна ковариации между s с тильдой и той же самой функцией от r.

То есть получается, что s с тильдой и s очень похожи, s с тильдой и s невозможно отличить, если смотреть только на математическое ожидание или на ковариацию с r, с r-квадрат, с любой функцией от r

Условное мат. ожидание — это Верно! $E(y_i | x_i)$ это наилучший прогноз $y_i$ формулируемый с помощью $x_i$.

На практике

Условия мат. ожидания

Пример

Здесь показано как рассчитывать

Условия мат. ожидания

Условное дисперсия

Если величины непрерывны и есть совместная функция плотности

Условия мат. ожидания

Свойства условного ожидания

Условия мат. ожидания

Найдите $E(x_i^2+2x_i | x_i)$ это $x^{2}_i+2x_i$. Потому что если мы знаем $x_i$, то и всё выражение слева мы легко спрогнозируем.

Условная дисперсия и условная ковариация

Условия мат. ожидания

Свойства:

Условия мат. ожидания

Пример

Упростите $Var(x^{2}_i+2x_i+\varepsilon_i|x_i)$ это равно $Var(\epsilon_i | x_i)$ Потому что Вспомним свойство $Var(s+h(r) | r)=Var(s|r)$ В роли s выступает $\varepsilon_i$. А если интуитивно: хотим спрогнозировать выражение слева, зная $x_i$ часть $x_i^2+2x_i$ нам доподлинно известна и на точность прогнозирования (а именно её меряет условная дисперсия) не влияет.

Геометрическая иллюстрация условного мат.ожидания

Запомнить

Интерпретируем дисперсию как квадрат длины случайной величины,
Интерпретируем корреляцию между случайными величинами как косинус угла между ними.

С помощью условного математического ожидания мы сформулируем стандартные предпосылки на случайную составляющую $\varepsilon$, а миеннто три предпосылки

Математическое ожидание от каждой случайной составляющей при известных иксах, то есть при известном каждом регресссоре для каждого наблюдения, я пишу коротко матрицу X. Условное математическое ожидание от $\varepsilon_i$ при всех известных регрессорах равна 0
Условное математическое ожидание от $\varepsilon^2$ при всех условных регрессорах равна $\sigma^2$. Или можно точно так же сказать, что условная дисперсия $\varepsilon_i$ при известной матрице X равна $\sigma^2$.
Kовариация между $\varepsilon_i$ и $\varepsilon__j$ при фиксированной матрице X равна нулю.

Предпосылки

Предпосылки про ковариацию и дисперсию можно коротко записать с помощью такого понятия как ковариационная матрица

Предпосылки

Когда мы говорим «ковариационная матрица некоего случайного вектора», мы имеем в виду здоровую табличку чисел. Первое число в первой строке — это дисперсия $\varepsilon_1$ . Второе число в первой строке, то есть первая строчка, второй столбец, (1,2) координаты, — это ковариация $\varepsilon_1$ и $\varepsilon_2$.

Соответственно, в ковариационной матрице, скажем, в третьей строчке, втором столбце находится ковариация $\varepsilon_3$ и $\varepsilon_2$, это третья строчка, второй столбец. С

оответственно, в этой матрице находятся и все дисперсии каждого $\varepsilon$ , и все попарные ковариации: $\varepsilon_1$ с $\varepsilon_3$ , $\varepsilon_2$ с $\varepsilon_7$ … все-все-все ковариации и дисперсии находятся в одной матрице, в одной табличке чисел.

Соответственно, первые наши две предпосылки на ε можно как сформулировать? Что ковариации равны нулю, а на диагонали дисперсии равны $\sigma^2$. То есть матрица принимает такой простой вид.

Предпосылки

Соответственно, в нашем случае наши предпосылки можно записать как ковариационная матрица вектора $\varepsilon$ при фиксированных регрессорах X равна $\sigma^2$ умножить на эту самую единичную матрицу, которая обозначается буковкой I, сокращение от английского «identity».

Итоги

У нас есть предпосылки, что дисперсия $\varepsilon$ при фиксированных регрессорах X равна $\sigma^2$ умножить на единичную матрицу, что на самом деле просто означает, что дисперсия $\varepsilon_i$ при фиксированном X равна $\sigma^2$, а ковариация разных $\varepsilon$ при фиксированном X равна нулю.
И у нас есть линейная модель,

Итоги

Это для примера двух объясняющих переменных. И, соответственно, эти предпосылки позволяют посчитать дисперсию любой оценки МНК $\beta_j$ с крышкой и любую ковариацию $\beta_j$ с крышкой и $\beta_l$ с крышкой.

Мы посчитаем для начала дисперсию, условную, оценки метода наименьших квадратов и ковариацию оценки метода наименьших квадратов для случая парной регрессии.

Итоги

Условная дисперсия МНК оценок. Доказательство

Условия нашей задачи

Что мы получаем. Запишем примечания

Примечания

Посчитаем в рамках предположений ковариацию $\widehat\beta_2$ и среднее

ПОлучаем не что иное как

Примечания

Итого в парной регрессии мы имеем

Итого в парной регрессии

Если регрессор $z$ сильно коррелирован с другими регрессорами, то Величина $RSS_z$ будет примерно равна 0 и поэтому дисперсия $Var(\widehat\beta_z|X)$ будет большой

Теорема

Штрих — это транспонированная матрица.

Свойства

Доказательство формулы для ковариационной матрицы

Этот фрагмент особо полезен тем, кто знает линейную алгебру. Оказывается, средствами линейной алгебры легко не только сразу посчитать дисперсию $\widehat\beta_1$ с крышкой или $\widehat\beta_2$ с крышкой, одним махом можно найти все эти дисперсии и все ковариации.

Свойства ковариационных матриц

Свойства

В этой матрице находятся и дисперсии условные каждого $\widehat\beta$, и ковариации каждого $\widehat\beta$ с каждым $\widehat\beta$. И в скалярном виде мы их выводили долго по отдельности и только для случая парной регрессии, а в матричном виде они выводятся легко, все и сразу, и этой формулой можно будет пользоваться.

Свойства

Оценка ковариационной матрицы

Константа $\sigma^2$ неизвестна, и именно эта константа входит в формулу для условной дисперсии любого оценённого коэффициента. А нам хочется строить доверительные интервалы для коэффициентов, проверять гипотезы, поэтому нам нужен какой-то способ оценить неизвестную константу.

Свойства

Статистические свойства оценок коэффициентов

Оценки коэффициентов, которые мы получаем методом наименьших квадратов, обладают рядом замечательных статистических свойств. Эти свойства, их очень много, и поэтому, чтобы как-то их было проще все осознать, мы их поделим на три части

Свойства

Предпосылки для применения метода наименьших квадратов и для получения хороших свойств коэффициентов.

Предпосылки

Свойства

Предположения на $\varepsilon_i$

Свойства

Предпоссылки на регрессоры

Свойства

Когда все эти предпосылки, которые мы сформулировали:

предпосылки о зависимости y от x,
предпосылки о распределении $\varepsilon$,
предпосылки о регрессорах,

когда все они выполнены, мы получаем, что верны три группы свойств.

Базовые

Свойства

Это очень хорошее свойство. Это говорит о том, что наш метод, конечно, может ошибаться, и оценка $\widehat\beta$, которую мы получаем, может не совпадать с настоящим $\beta$, но $\widehat\beta$ иногда будет больше настоящего $\beta$, иногда будет меньше настоящего $\beta$, но в среднем,$\widehat\beta$ попадает то влево, то вправо, но в среднем попадает в неизвестный коэффициент $\beta$

Свойства

Это очень хорошее свойство. Оно говорит о том, что если вы хотите простую оценку, то есть линейную, хотите оценку несмещенную, которая в среднем бы попадала в неизвестную истину, то ничего лучше оценок метода наименьших квадратов у нас не получится. То есть математически это означает, что условная дисперсия$\widehat\beta$ при фиксированных регрессорах альтернативного больше либо равна, чем условная дисперсия $\widehat\beta$, полученная по методу наименьших квадратов опять же при фиксированных регрессорах.

Свойства

То есть это свойство говорит, что оценка $\widehat\sigma^2$, предложенная нами, — RSS, делённая на (n – k), — она тоже несмещённая, и она несмещённым образом оценивает $\sigma^2$.

Асимптотические

Это те, что предполагают что число экспериментов велико

Свойства

Требующие нормальности

Свойства

Построение доверительных интервалов и проверка гипотез

Проверять гипотезы можно в двух случаях

Число наблюдений велико
Случайные ошибки нормальны

Построение доверительного интервала

Свойства

С ростом стандартной ошибки ширина доверительного интервала для коэффициента увеличивается

Проверка гипотез

Описание любого теста

Предпосылки теста
$H_0$ — проверяемая гипотеза
$H_1$ или $H_a$ — альтернативная или другими словами конкурирующая гипотеза
Формула для вычисления статистики
Закон распределения статистики при верной $H_0$

Последовательность действий

Формулируем гипотезу $H_0$ и выбираем уровень значимости. Это вероятность ошибки первого рода или вероятность отвергнуть $H_0$ при условии, что она верна. $\alpha = P (отвергнуть H_0| принять H_0)$. Можно формулировать гипотезу и под неё собирать данные
Рассчитываем наблюдаемое значение тестовой статистики $S_obs$
Находим критическое значение тестовой статистики $S_cr$
1. Сравниваем $S_obs$ и $S_cr$ и делаем вывод о $H_0$ — устаревший вариант. Сейчас всё чаще делают так:
1. Сравниваем P-значение и $\alpha$, делаем вывод о $H_0$. Если P-значение больше чем уровень $\alpha$, то гипотеза не отвергается. Если полученной P-значение меньше, то гипотеза отвергается

Пример. Доверительный интервал для коэффициента бета

считается это всё вот так

Пример

qt() — это функция в R

Пример. Доверительный интервал для дисперсии

Мы опираемся на теорему RSS делённое на сигма квадрат

qchisq() — это функция в R

Пример

Пример. Проверка гипотезы о коэффициенте beta

Предпологаем что коэффициент при доле сельскохозяйственного мужского населения равен нулю — другими словами фертильность не зависит от показателя того, насколько этот регион является сельскохозяйственным. Есть три способа проверить эту гипотезу

Пример

В литературе, скажем в научных статьях, очень часто стандартные ошибки выписывают под коэффициентами.

Интерпретация стандартной таблички

При оценке линейной модели регрессии, любой статистический пакет выдает более-менее стандартную табличку. Вот такую табличку выдает R.

Пример

Первый столбик в ней это, собственно, оценки коэффициентов. Т.е. смысл этого столбика, что мы можем записать уравнение линейной регрессии по используя эти коэффициенты.
Второй столбик — это стандартные ошибки. Это корни из диагональных элементов ковариационной матрицы
Третий столбик — это компьютер автоматом проверяет гипотезу о том, что на самом деле, зависимости от данной переменной нет. Он делает это при помощи Т-статистики. Т.е. это первый столбец делить на второй. T-статистика, проверяющая гипотезу о незначимости коэффициента может принимать любое значение. Знак T-статистики определяется знаком оценки коэффициента, и по модулю она может быть произвольной
P-value

Особенности проверки гипотез

Если $H_0$ не отвергается, это говорит о том, что зависимости нет.
Надо говорить аккуратно: «$H_0$ не отвергается» — это означает, что данные не противоречат гипотезе $H_0$
Значимость и существенность. Значимость это статистический факт — равен коэффициент нулю или не равен. Существенность — это на сколько он не равен нулю.
Стандартизировать коэффициенты перед анализом. Чтобы получить одинаково интерпритируемые единицы измерения
К сожалению, очень часто распространена такая порочная практика, что исследователь берет, включает кучу, не задумываясь о теоретической модели, включает кучу объясняющих переменных в свою модель и выбирает те из них, которые по t-статистикам оказались значимы. Это подход неправильный, поскольку как только мы согласились на некую вероятность ошибки первого рода, например, мы выбрали вероятность ошибки первого рода типичную в экономических приложениях 5%. «Я запустил регрессию на кучу переменных и отобрал те, которые значимы» — это неправильный подход.

Проверка гипотезы о связи коэффициентов

К примеру, мы хотим в рамках нашей модели проверить, что воздействие, рост фертильности, вызванный увеличением доли мужчин, занятых в сельском хозяйстве, одинаков по силе с ростом фертильности при росте католического населения, то есть я хочу проверить гипотезу о том, что разница этих двух коэффициентов равна нулю. Как это сделать?

Гипотеза не отвергается

Второй способ — подобрать коэффициенты. Т.е. отнять и прибавить коэффициент при втором регрессоре. Программа автоматом рассчитает вероятность того, что разница коэффициентов равна нулю

Гипотеза не отвергается

Задания на R

Подключаем нужные пакеты

library(memisc)

## Loading required package: lattice

## Loading required package: MASS

## 
## Attaching package: 'memisc'

## The following object is masked from 'package:ggplot2':
## 
##     syms

## The following objects are masked from 'package:stats':
## 
##     contr.sum, contr.treatment, contrasts

## The following object is masked from 'package:base':
## 
##     as.array

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

## ✔ tibble  2.1.3     ✔ purrr   0.3.2
## ✔ tidyr   0.8.2     ✔ dplyr   0.8.3
## ✔ readr   1.3.1     ✔ stringr 1.4.0
## ✔ tibble  2.1.3     ✔ forcats 0.4.0

## ── Conflicts ────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::%@%()     masks memisc::%@%()
## ✖ dplyr::collect() masks memisc::collect()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ dplyr::recode()  masks memisc::recode()
## ✖ dplyr::rename()  masks memisc::rename()
## ✖ dplyr::select()  masks MASS::select()
## ✖ dplyr::syms()    masks memisc::syms(), ggplot2::syms()

library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(sjPlot)
library(sgof)

## Loading required package: poibin

library(foreign)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:memisc':
## 
##     recode

library(hexbin)
library(rlms)

Работа со случайными величинами

Генерация случайных величин.

100 величин, распределение нормально. Любое распределение генерится в R функцией, которая начинается с r и далее название распределения.

z <- rnorm(100, mean = 5, sd = 3)
z[56]

## [1] 4.414769

z[2:9]

## [1] 6.986693 3.969756 1.714242 3.768766 3.898127 8.314513 6.326293 8.544909

qplot(z)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Построим функцию плотности

Все функции плотности начинаются с буквы d (от density) в R.

x <- seq(-10, 15, by = 0.5)
y <- dnorm(x, mean = 5, sd = 3)
qplot(x, y, geom = "line")

Расчёт вероятностей

Все функции плотности начинаются с буквы p (от probability) в R.

pnorm(3, mean = 5, sd = 3)

## [1] 0.2524925

если я хочу найти вероятность того, что z лежит в диапазоне от 4 до 9, то это есть с точки зрения здравого смысла вероятность того, что z меньше 9 минус вероятность того, что z меньше 4

pnorm(9, mean = 5, sd = 3) - pnorm(4, mean = 5, sd = 3)

## [1] 0.5393474

Квантили распределения. $P(Z<a)=0.7$ какое $a$

Все функции плотности начинаются с буквы q (от quantile) в R.

qnorm(0.7, mean = 5, sd = 3)

## [1] 6.573202

Квантиль — это на самом деле, обратная функция к функции распределения. То есть, если, например, я хочу найти такое число а, чтобы вероятность того, что z меньше а была равна, скажем, 0.7. И вот надо найти такое а. То соответственно, это можно найти с помощью квантильной функции

Есть разные популярные распределения. chisq — хи квадрат. t — стьюдента. f — Эф распределение

Проверка гипотез о распределениях

Множественная регрессия, проверка гипиотез

h <- swiss
glimpse(h)

## Observations: 47
## Variables: 6
## $ Fertility        <dbl> 80.2, 83.1, 92.5, 85.8, 76.9, 76.1, 83.8, 92.4,…
## $ Agriculture      <dbl> 17.0, 45.1, 39.7, 36.5, 43.5, 35.3, 70.2, 67.8,…
## $ Examination      <int> 15, 6, 5, 12, 17, 9, 16, 14, 12, 16, 14, 21, 14…
## $ Education        <int> 12, 9, 5, 7, 15, 7, 7, 8, 7, 13, 6, 12, 7, 12, …
## $ Catholic         <dbl> 9.96, 84.84, 93.40, 33.77, 5.16, 90.57, 92.85, …
## $ Infant.Mortality <dbl> 22.2, 22.2, 20.2, 20.3, 20.6, 26.6, 23.6, 24.9,…

Мы оценим линейную модель регрессии. Будем предполагать, что фертильность, Fertility, зависит от доли католического населения в данном кантоне, от показателя, насколько это регион сельскохозяйственный, и, скажем, от Examination.

model <- lm(data = h, Fertility ~ Catholic + Agriculture + Examination)

Посмотрим отчет по этой моделе

summary(model)

## 
## Call:
## lm(formula = Fertility ~ Catholic + Agriculture + Examination, 
##     data = h)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.746  -5.724   1.248   6.324  18.849 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 90.86803    8.63691  10.521 1.81e-13 ***
## Catholic     0.04240    0.04148   1.022 0.312347    
## Agriculture -0.09516    0.08587  -1.108 0.273930    
## Examination -1.07035    0.27316  -3.918 0.000315 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.616 on 43 degrees of freedom
## Multiple R-squared:  0.4461, Adjusted R-squared:  0.4074 
## F-statistic: 11.54 on 3 and 43 DF,  p-value: 1.113e-05

Отсюда можем сказать, что гипотеза о том, что $\beta_1$ равно нулю, отвергается; гипотеза о том, что $\beta_2$ равно нулю, не отвергается; гипотеза о том, что $\beta_3$ равно нулю, не отвергается; гипотеза о том, что $\beta_4$ равно нулю, отвергается.

Только значения коэффициентов, И также можем легко получить доверительные интервалы:

coeftest(model)

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) 90.868026   8.636909 10.5209 1.807e-13 ***
## Catholic     0.042401   0.041476  1.0223 0.3123467    
## Agriculture -0.095159   0.085867 -1.1082 0.2739301    
## Examination -1.070347   0.273160 -3.9184 0.0003147 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

confint(model)

##                   2.5 %       97.5 %
## (Intercept) 73.45003853 108.28601395
## Catholic    -0.04124219   0.12604494
## Agriculture -0.26832536   0.07800811
## Examination -1.62122538  -0.51946784

Проверка линейных гипотез

проверим линейную гипотезу о том, что коэффициент зависимости при доле католического населения и при доле населения, занятого в сельском хозяйстве, одинаковые.

Воспользуемся хитрым способом, в котором мы складываем две переменные — способ с построением вспомогательной регрессии — чтобы проверить гипотезу.

Значок I означает инструкцию для R, что надо трактовать Catholic + Agriculture — «плюс» в прямом смысле.

model_aux <- lm(data = h, 
                Fertility ~ Catholic + I(Catholic + Agriculture) + Examination)

Выведем отчёт о модели

summary(model_aux)

## 
## Call:
## lm(formula = Fertility ~ Catholic + I(Catholic + Agriculture) + 
##     Examination, data = h)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.746  -5.724   1.248   6.324  18.849 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               90.86803    8.63691  10.521 1.81e-13 ***
## Catholic                   0.13756    0.09585   1.435 0.158483    
## I(Catholic + Agriculture) -0.09516    0.08587  -1.108 0.273930    
## Examination               -1.07035    0.27316  -3.918 0.000315 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.616 on 43 degrees of freedom
## Multiple R-squared:  0.4461, Adjusted R-squared:  0.4074 
## F-statistic: 11.54 on 3 and 43 DF,  p-value: 1.113e-05

Коэффициент при Catholic незначим, потому что P-уровень равен 0.158483. Это говорит о том, что гипотеза о том, что коэффициенты при Catholic и Agriculture равны, не отвергается.

Таким образом, мы смогли проверить гипотезу о том, что два коэффициента истинных, неизвестных нам, равны, и эта гипотеза в нашем случае не отвергается.

Стандартизированные коэффициенты и эксперимент с ложно-значимыми регрессорами

Один из способов почувствовать существенный коэффициент или нет, — это посчитать стандартизированные коэффициенты $\widehat\beta$, то есть привести все объясняющие переменные и объясняемую переменную к неким универсальным единицам измерения, чтобы они были сравнимы, а именно: вычесть из каждой переменной её среднее и поделить на оценённое стандартное отклонение.

Шаг 1 — Стандартизация коэффициентов

На шаге один мы преобразуем каждую переменную, масштабируем каждую переменную. Значит, создадим набор h_st, где мы изменим каждую переменную, функция называется mutate_each, в наборе данных h, а формула, по которой мы будем, функция, по которой мы будем менять каждую переменную, называется scale, она осуществляет как раз масштабирование, вычитание среднего и деления на стандартную ошибку.

h_st <- mutate_each(h, "scale")
glimpse(h_st)

## Observations: 47
## Variables: 6
## $ Fertility        <dbl[,1]> <matrix[25 x 1]>
## $ Agriculture      <dbl[,1]> <matrix[25 x 1]>
## $ Examination      <dbl[,1]> <matrix[25 x 1]>
## $ Education        <dbl[,1]> <matrix[25 x 1]>
## $ Catholic         <dbl[,1]> <matrix[25 x 1]>
## $ Infant.Mortality <dbl[,1]> <matrix[25 x 1]>

Шаг 2 — построение модели

model_st <- lm(data = h_st, Fertility ~ Catholic + Agriculture + Examination)
summary(model_st)

## 
## Call:
## lm(formula = Fertility ~ Catholic + Agriculture + Examination, 
##     data = h_st)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.14113 -0.45821  0.09989  0.50623  1.50891 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.729e-16  1.123e-01   0.000 1.000000    
## Catholic     1.416e-01  1.385e-01   1.022 0.312347    
## Agriculture -1.730e-01  1.561e-01  -1.108 0.273930    
## Examination -6.836e-01  1.745e-01  -3.918 0.000315 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7698 on 43 degrees of freedom
## Multiple R-squared:  0.4461, Adjusted R-squared:  0.4074 
## F-statistic: 11.54 on 3 and 43 DF,  p-value: 1.113e-05

Шаг 3 — визуализация коэффициентов.

Шаг 4 — принятие решения

Искусственный эксперимент

Сейчас мы на искусственных данных проиллюстрируем идею того, что нельзя просто так построить регрессию и отвечать на вопрос, а какие коэффициенты у меня значимы.

Мы сочиним в нашем искусственном эксперименте совершенно несвязанный игрек, который никак не зависит от якобы объясняющих переменных. У нас будет 40 якобы объясняющих переменных, одна якобы зависимая, хотя на самом деле независимая, и мы будем строить, оценивать модель линейной регрессии.

Генерим данные

set.seed(42)
d <- matrix(rnorm(100*41, mean = 0, sd = 1), nrow = 100)
df <- data.frame(d)
glimpse(df)

## Observations: 100
## Variables: 41
## $ X1  <dbl> 1.37095845, -0.56469817, 0.36312841, 0.63286260, 0.40426832,…
## $ X2  <dbl> 1.200965376, 1.044751087, -1.003208647, 1.848481902, -0.6667…
## $ X3  <dbl> -2.00092924, 0.33377720, 1.17132513, 2.05953924, -1.37686160…
## $ X4  <dbl> -0.0046207678, 0.7602421677, 0.0389909129, 0.7350721416, -0.…
## $ X5  <dbl> 1.33491259, -0.86927176, 0.05548695, 0.04906691, -0.57835573…
## $ X6  <dbl> 1.029140719, 0.914774868, -0.002456267, 0.136009552, -0.7201…
## $ X7  <dbl> -0.2484829, 0.4223204, 0.9876533, 0.8355682, -0.6605219, 1.5…
## $ X8  <dbl> 0.29469236, 0.39274127, -1.00084371, -0.32572712, -1.0083488…
## $ X9  <dbl> 0.68880775, 0.72508302, 0.21738021, -0.20165673, -1.36568986…
## $ X10 <dbl> 0.94192422, -0.24861404, 0.09647886, -0.43393094, 2.17866787…
## $ X11 <dbl> 2.32505849, 0.52412218, 0.97073342, 0.37697340, -0.99593340,…
## $ X12 <dbl> -0.6502844, -1.0031832, -0.5351139, -0.1104146, 0.6004302, 0…
## $ X13 <dbl> -0.74651645, 0.03660612, 0.32330962, 0.37967603, 0.87655650,…
## $ X14 <dbl> -1.18554694, -0.65838929, 1.08950795, 0.50878594, -0.1359066…
## $ X15 <dbl> 0.87729465, -1.77337144, -0.04568732, -0.39487225, -0.128056…
## $ X16 <dbl> -0.60138300, -0.13581614, -0.98727278, 0.83192501, -0.795059…
## $ X17 <dbl> 0.02931665, 1.93431434, 0.83751088, -0.18534633, 0.53332918,…
## $ X18 <dbl> 1.44820167, -0.08585566, -1.57376980, 1.29247616, 0.08446776…
## $ X19 <dbl> -1.2213334, -0.4529860, -0.7023199, -0.6743804, -2.3030232, …
## $ X20 <dbl> 0.09782704, 1.46502103, 0.99043585, 1.41438675, -1.33589777,…
## $ X21 <dbl> 0.25057807, -0.27792405, -1.72473573, -2.00670494, -1.291808…
## $ X22 <dbl> -0.58693858, 0.92697573, -0.06540585, -0.93711176, -0.633859…
## $ X23 <dbl> 0.19033984, -0.07173877, -0.00285171, -1.10821896, 0.9351917…
## $ X24 <dbl> 0.75806083, 1.47961434, -0.61956241, 0.25968723, -1.08652092…
## $ X25 <dbl> 1.32750456, -0.60083583, 0.05650686, -0.53107629, -0.0808987…
## $ X26 <dbl> 0.617336710, -0.004541141, -0.091256360, 0.399959375, 0.5889…
## $ X27 <dbl> -0.2418159, 0.5721363, 0.7057918, 0.2048996, 0.5521034, -0.5…
## $ X28 <dbl> -0.54468934, -0.43453296, -0.14084519, 1.75814127, -0.830141…
## $ X29 <dbl> -0.8495311, 0.1151387, -1.3313372, 2.3744583, -2.4958364, 2.…
## $ X30 <dbl> -0.63258679, 1.06719020, 0.33739592, -0.95364566, 0.25117357…
## $ X31 <dbl> -0.68566166, -0.79271447, -0.40700415, -1.14867061, 1.115760…
## $ X32 <dbl> 1.37843190, -2.95313731, 0.95583656, 1.06003423, 1.35096057,…
## $ X33 <dbl> -2.0645260, 0.6772492, 0.8363831, 0.1222439, -0.7882496, 0.6…
## $ X34 <dbl> -0.9762936, -1.0083745, -0.2563995, 0.8828167, 0.3760271, 0.…
## $ X35 <dbl> -0.1062227, -1.2870865, -0.9080599, -0.2839357, 0.2665925, 1…
## $ X36 <dbl> 0.29043335, -1.30379423, -0.11410060, -0.09296209, -0.122560…
## $ X37 <dbl> -0.18434938, -0.26304466, -0.52564624, -0.67522750, -0.15678…
## $ X38 <dbl> -0.59179607, -0.88500175, -0.52743656, 0.03103741, -0.946515…
## $ X39 <dbl> -0.007239089, -0.462310513, -0.724387536, -0.639016920, 0.44…
## $ X40 <dbl> -0.4503987, -1.3772278, 1.0485915, 0.3055562, -0.9187486, -0…
## $ X41 <dbl> -0.1418087, -0.8138981, -0.3255406, 0.3781574, -1.9944854, -…

Будем объяснять первую переменную, всеми оставшимися в наборе данных переменными. Для этого есть удобное сокращение в виде точки .

model_pusto <- lm(data = df, X1~.)

summary(model_pusto)

## 
## Call:
## lm(formula = X1 ~ ., data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7476 -0.6018  0.1303  0.5730  1.5573 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.001317   0.143255   0.009   0.9927  
## X2           0.142212   0.161744   0.879   0.3828  
## X3          -0.314385   0.136043  -2.311   0.0244 *
## X4           0.031564   0.166652   0.189   0.8504  
## X5           0.178323   0.150741   1.183   0.2416  
## X6           0.111487   0.137823   0.809   0.4218  
## X7          -0.108252   0.137431  -0.788   0.4340  
## X8           0.145668   0.149928   0.972   0.3352  
## X9          -0.036021   0.143237  -0.251   0.8023  
## X10          0.110686   0.121645   0.910   0.3666  
## X11         -0.038508   0.144271  -0.267   0.7905  
## X12          0.034769   0.145851   0.238   0.8124  
## X13          0.049680   0.152415   0.326   0.7456  
## X14         -0.158988   0.142948  -1.112   0.2706  
## X15         -0.014374   0.148596  -0.097   0.9233  
## X16          0.248442   0.172522   1.440   0.1551  
## X17         -0.061389   0.144079  -0.426   0.6716  
## X18         -0.047658   0.138050  -0.345   0.7312  
## X19          0.235795   0.141078   1.671   0.0999 .
## X20         -0.034408   0.148917  -0.231   0.8181  
## X21          0.001259   0.131329   0.010   0.9924  
## X22         -0.180049   0.154533  -1.165   0.2487  
## X23         -0.072289   0.132548  -0.545   0.5875  
## X24         -0.032262   0.190246  -0.170   0.8659  
## X25          0.052217   0.148228   0.352   0.7259  
## X26         -0.206708   0.139283  -1.484   0.1431  
## X27         -0.082317   0.121478  -0.678   0.5007  
## X28         -0.096665   0.164121  -0.589   0.5581  
## X29          0.202307   0.133146   1.519   0.1340  
## X30         -0.047584   0.129551  -0.367   0.7147  
## X31          0.055963   0.162809   0.344   0.7323  
## X32          0.307859   0.158290   1.945   0.0566 .
## X33          0.044117   0.150258   0.294   0.7701  
## X34         -0.093703   0.135581  -0.691   0.4922  
## X35         -0.344475   0.168765  -2.041   0.0457 *
## X36         -0.065275   0.135642  -0.481   0.6321  
## X37          0.063727   0.143506   0.444   0.6586  
## X38          0.132398   0.161363   0.820   0.4152  
## X39         -0.055207   0.158315  -0.349   0.7285  
## X40         -0.059821   0.174217  -0.343   0.7325  
## X41         -0.085292   0.142895  -0.597   0.5529  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.123 on 59 degrees of freedom
## Multiple R-squared:  0.3074, Adjusted R-squared:  -0.1622 
## F-statistic: 0.6546 on 40 and 59 DF,  p-value: 0.9211

У нас получилось, что переменные X3 и X35 ылияют на первую переменную при 5% уровне значимости, а X19 и X32 влияют на 10% уровне значимости.

С чем это связано? Это связано с тем, что когда мы фиксируем уровень значимости, мы соглашаемся на некоторую вероятность ошибиться.

Соответственно, когда мы зафиксировали уровень значимости 10%, это означает, что с вероятностью 10% мы в случае на самом деле отсутствия зависимости якобы её обнаружим.

Вывод

Соответственно, из этого искусственного эксперимента нужно сделать простой вывод, что стратегия «я оценю модель с большим количеством объясняющих переменных и выпишу из них значимые, и скажу, что от них игрек зависит», — неправильная, потому что в силу того, что есть для каждого регрессора десятипроцентный или пятипроцентный шанс сделать ошибку, при большом количестве регрессоров кто-то якобы значимым будет, даже если на самом деле никакой зависимости нет.

Сравнить несколько моделей

Есть несколько моделей

model <- lm(data = h, Fertility ~ Catholic + Agriculture + Examination)
model2 <- lm(data = h, Fertility ~ Catholic + Agriculture)

Как нам лучше их сравнить. Сделаем табличку, сравниваем при помощи функции mtable() из пакета memisc

compare_12 <- mtable(model, model2)

RLMS

beta <- c(500, 30)

Untitled

Неделя один

Обозначения

Пример

План действий

Метод наименьших квадратов

Случай для одного регрессора

Случай для двух регрессоров

Случай для множества регрессоров

Графическое представление

Множественные регрессоры

Суммы квадратов

Лекбез по линейной алгебре

Обозначения

Линейная модель регрессии в n-мерном пространстве

Геометрия множественной регрессии

Коэффициент детерминации

Пример с фертильностью. МНК в R

Неделя два

Условное математическое ожидание

Пример

Условное дисперсия

Свойства условного ожидания

Условная дисперсия и условная ковариация

Геометрическая иллюстрация условного мат.ожидания

Итоги

Условная дисперсия МНК оценок. Доказательство

Доказательство формулы для ковариационной матрицы

Оценка ковариационной матрицы

Статистические свойства оценок коэффициентов

Базовые

Асимптотические

Требующие нормальности

Построение доверительных интервалов и проверка гипотез

Построение доверительного интервала

Проверка гипотез

Описание любого теста

Последовательность действий

Пример. Доверительный интервал для коэффициента бета

Пример. Доверительный интервал для дисперсии

Пример. Проверка гипотезы о коэффициенте beta

Интерпретация стандартной таблички

Особенности проверки гипотез

Проверка гипотезы о связи коэффициентов

Задания на R

Работа со случайными величинами

Проверка гипотез о распределениях

Стандартизированные коэффициенты и эксперимент с ложно-значимыми регрессорами

Искусственный эксперимент

Вывод

Сравнить несколько моделей

RLMS